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Abstract 


Granular material is showing very often in geotechnical engineering, petroleum 
engineering, material science and physics. The packings of the granular material play 
a very important role in their mechanical behaviors, such as stress-strain response, 
stability, permeability and so on. Although packing is such an important research topic 
that its generation has been attracted lots of attentions for a long time in theoretical, 
experimental, and numerical aspects, packing of granular material is still a difficult 
and active research topic, especially the generation of random packing of non-spherical 
particles. To this end, we will generate packings of same particles with same shapes, 
numbers, and same size distribution using geometry method and dynamic method, 
separately. Specifically, we will extend one of Monte Carlo models for spheres to 
ellipsoids and poly-ellipsoids. 
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1 Introduction 


Granular material is showing very often in geotechnical engineering, petroleum engineering, 
material science and physics [Buchalter and Bradley, 1994, Jaeger et al., 1996, Regueiro 
et al., 2014b, Senseney et al., 2017, Zhang, 2016]. The packings of the granular material play 
a very important role in their mechanical behaviors, such as stress-strain response, stability, 
permeability and so on [Troadec et al., 1991, Maggi et al., 2008, Zhang and Regueiro, 2015, 
Amirrahmat et al., 2018, Regueiro et al., 2014a, Luo et al., 2022]. For example, Troadec 
et al. [1991] studied experimentally the effects of different packings on the stress-strain 
relation during compression of packed cylinders and found that the coordination number of 
the packing govern the stress-strain paths during compression. Bassett et al. [2012], Sadd 
et al. [2000] then further concluded experimentally and numerically, respectively, that force 
chain of packed granular materials plays significant effects on sound propagation in granular 
materials. And our previous researches [Zhang et al., 2015] show that different packed sand 
grains behave differently in the simulations of Split Hopkinson Pressure Bar experiments. 
Although packing is such an important research topic that its generation has been at- 
tracted lots of attentions for a long time in theoretical [Furnas, 1929, Melissen and Schuur, 
1995, Brouwers, 2006, Torquato and Jiao, 2009], experimental [Furnas, 1929, Bideau and 
Hansen, 1993], and numerical [He et al., 1999, Jia and Williams, 2001, Jiang et al., 2003, 
Dong et al., 2006, Zhou et al., 2011] aspects, packing of granular material is still a difficult 
and active research topic, especially the generation of random packing of non-spherical par- 
ticles. There are two kinds of numerical generation methods of random packing based on 
geometry and dynamic, respectively. Geometry method is moving and rotating particles in 
a packing by geometry constraints to reduce overlap between particles. After overlap is de- 
creased to a small tolerance, the packing is viewed as a stable system and packing generation 
is finished. The geometry method is efficient and has been applied successfully to generate 
random packings of monodisperse spheres, polydisperse spheres, ellipsoids, super-ellipsoids, 


and spherocylinders [Buchalter and Bradley, 1994, Yang et al., 1996, He et al., 1999, Jia and 


Williams, 2001, Kansal et al., 2002, Fu and Dekelbab, 2003, Abreu et al., 2003, Donev et al., 
2007, Maggi et al., 2008, Stafford and Jackson, 2010, Delaney and Cleary, 2010]. Among 
so many geometry packing generation algorithms, Monte Carlo model is very famous and 
extended to many different variables. The dynamic method is simulating the physical inter- 
actions between particles during packing process [Cheng et al., 2000, Jiang et al., 2003, Zhou 
et al., 2011] based on Discrete Element Method (DEM). DEM [Cundall and Strack, 1979] 
is used to simulate the motion of individual particles/grains within an overall deformation 
and/or flow of a granular medium, which is borrowed to simulate the dynamic process during 
packing generation. For example, Cheng et al. [2000] randomly generated circular particles 
inside a box without any overlaps by removing the overlapped particles, then these random 
particles are settled down by gravity force to generate random packings. Jiang et al. [2003] 
generated particles randomly in a container without any overlaps, and then compress these 
particles to reach specific packing fraction. Zhou et al. [2011] generated packings by pouring 
ellipsoidal particles from a certain height into a container, which is more physically related 
to the generation of packings in experiments. These descriptions of the two methods suggest 
that geometry methods are much more faster than dynamic methods, since it takes a lot of 
time for the particles in the latter to come to rest. While dynamic methods considered the 
interactions between particles mechanically, and the pouring method is more similar phys- 
ically to the generation of packings in experiments. Although both geometry methods and 
dynamic methods can generate packings with similar fractions as in experiments, there are 
worries about what it will influence on the packings if geometry methods do not consider 
packing as a dynamic process involving forces between particles [Jia and Williams, 2001, 
Dong et al., 2006, An et al., 2008, Zhou et al., 2011]. The motivation of this paper is then to 
explore the influence and to answer if we can take advantages of the efficiency of geometry 
methods without introducing much difference to the final packings of particles compared 
with dynamic generating methods. To this end, we will generate packings of same particles 


with same shapes, numbers, and same size distribution using geometry method and dynamic 


method, separately. Specifically, we will extend one of Monte Carlo models proposed by He 
et al. [1999] for spheres to ellipsoids and poly-ellipsoids [Peters et al., 2009]. To simulate the 
dynamic process during packings, we will adopt gravitational deposition [Yan et al., 2010] 
of particle assembly which is similar as the pouring method described in [Zhou et al., 2011]. 
After packings are generated by the two methods separately, the corresponding packings 
will be analyzed and compared with respect to packing metrics, such as packing fraction, 
Coordination Number (CN), Radia Distribution Function (RDF) [Konakawa and Ishizaki, 
1990, He et al., 1999, Zhou et al., 2011], fabric tensor [Yimsiri and Soga, 2010, 2011], and 
also force chain [Sadd et al., 2000, Bassett et al., 2012]. 

Currently, most of packings are generated for spheres and ellipsoids. Some researchers 
studied and concluded that particle shape plays a crucial role in packings of ellipsoids 
[Bezrukov and Stoyan, 2006], super-ellipsoids [Delaney and Cleary, 2010], hyperspheres 
[Skoge et al., 2006], and sphereocylinders [Abreu et al., 2003]. Thus this paper also consid- 
ered packings of different shapes generated by the two methods, such as spheres, ellipsoids 
(prolate and oblate) [Bezrukov and Stoyan, 2006], and poly-ellipsoids (carrot and half-dome) 
[Peters et al., 2009]. With poly-ellipsoids, we can create many unsymmetrically shaped par- 
ticles, it is very important since most of geotechnical material are not regular or symmetric, 


and no one else has studied packings of unsymmetric particles. 


2 Algorithms 


This paper considered two packing generation methods, an extended Monte Carlo Model and 
DEM simulation of gravitational deposition, which are belonging to the two major packing 
generation methods in community, respectively. The two methods will be used to generate 
packings of same particle assemblies, and then the packings will be analyzed and compared 
to study the effects caused by not considering dynamic interactions between particles in 


Monte Carlo packing algorithm. 


2.1 Extended Monte Carlo Model 


He et al. [1999] presented a variable of Monte Carlo model to generate random packings of 
unequal spherical particles obeying any specified distributions. He et al. [1999] first generated 
a certain number (n) of spherical particles following a specified size distribution, and these 
particles are uniformly randomly placed within a cubic domain with the initial size as Lo, 


which is determined as [He et al., 1999] 


pa q” 
De g (1) 


where ®p is an arbitrary initial packing density which can be unphysically higher than 
random close packing fraction, in this paper ®o is taken as 0.86 following He et al. [1999]’s 
suggestion, and v; is the volume of particle p;. The meaning of Eqn.1 is generating a cubic 
domain containing the initial particle assembly with a packing fraction as ®o. 

There will be lots of overlaps between particles in this initial packing. During each step, 
He et al. [1999] moved an overlapped particle by the sum vector considering all overlaps of 
this particle. For example, the new position of particle p; with n; overlapped particles is 


given by [He et al., 1999] 
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where t;; is the translation vector of particle p; caused by the overlapping with particle p,, 
which is shown in Fig.1. After a certain number of iterations, the mean of the relative overlap 
will eventually drop to a very small preset tolerance, such that the packing is regarded as a 
stable packing. 

To extend this Monte Carlo model for ellipsoids and poly-ellipsoids, the major changes are 
the contact detection and rotations of overlapped particles. The contact detection algorithm 


is to find contacts between particles and their overlaps. For spheres, the contact detection 


Figure 1. Push outward particle p; in the direction of branch vector connecting the centroids of 
particle p; and p; by translation vector t;; in order to eliminate the overlap between particle p; and 
pj, motivated by He et al. [1999]. 

algorithm is trivial, however it is not straightforward for ellipsoids and poly-ellipsoids, the 
contact detection algorithms for which are described in Yan et al. [2010] and Zhang et al. 
[in review], respectively. The other change is to add rotations of overlapped particles for 
ellipsoids and poly-ellipsoids. Suppose we have two overlapping particles (p; and p;) as in 
Fig.2, then the translation and rotation of particle p; caused by particle p; in our model are, 


respectively 


tij = TA (3) 
where m;, mj are the mass for particles p; and pj, and A, d; are shown in Fig.2. A is the 


overlap vector pointing from point 7 to point 7, which are contact points on the surface of 


particle p; and pj, respectively. d; is a vector pointing from centroid p; to point 7. Similarly, 


if particle p; has n; overlapped particles, then the overall translation vector and rotation 


vector are 


(5) 


(6) 


Then after a certain number of steps, the mean of relative overlap will drop below a tolerance 


and the packing generation is finished. In this research, 1 x 1074 is accepted as the tolerance 


for the mean of relative overlap. 


Figure 2. 2D illustration of overlapped particles. A is the overlap vector connecting two contact 
points 2 and j, which are on the surface of particle p; and pj. d; is pointing from centroid p; to 


point 7. 


2.2 DEM simulation of gravitational deposition 


DEM [{Cundall and Strack, 1979] is often used in geotechnical engineering to simulate the 


motion of individual particles/grains within a deforming/flowing granular medium|Horner 


et al., 2001, Hopkins, 2004, Zhu et al., 2008, Yan et al., 2010, O’Sullivan, 2011, Knuth et al., 


2012]. According to Cundall and Strack [1979], Yan et al. [2010], the governing equation for 


the translation and rotation of DEM particle 7 in an assembly are 


m,u; = F; 
, (7) 

L0; = M; 
where u is the particle displacement; 0, the orientational vector of the particle; m, the particle 
mass; J, the moment of inertia of the particle; F, the resultant force, which includes body 
force and contact forces by its overlapping particles; and M, the resultant moment about the 
principal axes of inertial frame. In this DEM model, Hertz-Mindlin contact theory [Mindlin, 
1949] is applied which includes nonlinear elasticity and slip, with addition of Coulomb friction 
model to evaluate tangential stick-slip conditions. Central difference time integration method 

is used to solve Eqn.7. 

To generate packings of granular material, a gravitational deposition of particle assembly 
is simulated dynamically using DEM. The particles will be positioned inside a container with 
open top boundary. The particles are positioned initially without any overlaps with other 
particles and boundaries. Then these particles will be dropped by the gravity forces at zero 


initial velocities. After all the particles come to rest, the packing generation is completed. 


2.3 Initial particle size distribution 


Similar as He et al. [1999], sizes of particles obey truncated log-normal distribution, the 


probability density function of particle with radius r is 


f(r) = —— e- tnr-tnro)?/ (20) (8) 

y 2r)or 
In this research, the mean radius is chosen as rọ = 1m, and the standard deviation ¢ = 
0.25m, and we only allow particle radius between 0.2m and 2.5m. The particle size distri- 


bution generated in numerical sample is compared with analytical log-normal distribution 


as shown in Fig.3. 
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Figure 3. Particle size distribution in numerical sample is compared well with the analytical 
log-normal distribution. 


There are five types of shaped particles in this research, i.e. sphere, prolate, oblate, 
carrot, and half-dome. Prolate and oblate are ellipsoids, while carrot and half-dome are 
poly-ellipsoids. For spheres, their radii are determined by Eqn.8. For ellipsoids and poly- 
ellipsoids, we make their major semi-lengths equal to r determined by Egqn.8. Particle 
assembly with each shaped particle will be only generated once and then the same assembly 


will be used in the extended Monte Carlo model and DEM gravitational deposition. 
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